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Biochemical reaction networks in living cells usually involve reversible covalent modification of 
signaling molecules, such as protein phosphorylation. Under conditions of small molecule numbers, 
as is frequently the case in living cells, mass action theory fails to describe the dynamics of such 
systems. Instead, the biochemical reactions must be treated as stochastic processes that intrinsically 
generate concentration fluctuations of the chemicals. We investigate the stochastic reaction kinetics 
of covalent modification cycles (CMCs) by analytical modeling and numerically exact Monte-Carlo 
simulation of the temporally fluctuating concentration. Depending on the parameter regime, we find 
for the probability density of the concentration qualitatively distinct classes of distribution functions, 
including power law distributions with a fractional and tunable exponent. These findings challenge 
the traditional view of biochemical control networks as deterministic computational systems and 
suggest that CMCs in cells can function as versatile and tunable noise generators. 

PACS numbers: 87.18.Vf, 87.10.Mn, 82.20.Fd, 05.10.Gg, 82.20.Db, 87.15. R-, 82.20.-w, 82.37.-j, 82.39.-k, 
82.40.-g, 05.40.-a 



I. INTRODUCTION 

Living cells transduce chemical signals from the envi- 
ronment via trans-membrane receptors to their interior. 
The activated receptors trigger chains of chemical reac- 
tions along so-called signaling pathways, which can for 
example lead to the expression of selected genes in re- 
sponse to the external stimulus. Complex reaction net- 
works arise when several linear pathways are cross-linked 
by multiple biochemical interactions. Such signal trans- 
duction networks are traditionally thought of as deter- 
ministic " computers" , in which information is coded by 
the relative concentration of bio-chemicals. This study 
challenges this view and suggests that stochastic concen- 
tration fluctuations are the primary mode of operation 
for most intracellular signaling cascades. 

It is well known that the numbers of receptors and sig- 
naling molecules fluctuate as a function of time and from 
cell to cell |Fur05] . The role of these fluctuations, often 
regarded as noise, is still poorly understood. How can 
cells properly react to external stimuli when the signals 
have to pass through noisy channels ? Is the degree of 
noise actively suppressed for certain key signaling pro- 
teins ? Or is the present understanding of intra-cellular 
control, based on mass action theory, overly simplified ? 

Using covalent modification cycles (CMCs) as a sim- 
ple model system, we show that the magnitude of con- 
centration fluctuations, relative to the mean value, can 
indeed be enormous. We demonstrate that CMCs can 
be viewed as versatile and tunable noise generators. De- 
pending on the system parameters, qualitatively different 
classes of probability density functions (PDFs) of concen- 
tration fluctuations emerge, including extremely broad 
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and asymmetric distributions with fractional power law 
tails. 

CMCs are a very common motif in cellular reaction 
networ ks |Kit02l IKho06l IHar99l IKos981 ISha84l IKre811 
ISta77| . The typical structure of a CMC is shown in Fig[l] 
below. In such systems, a substrate protein is found in 
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FIG. 1: (Color online) Schematic of a covalent modification 
cycle (CMC): Substrate Xo is activated by enzyme A into 
the modified form X and deactivated by enzyme D. Each 
shaded sub-module denotes an enzymatic conversion reaction 
(Unbinding reactions not shown). 

two different chemical states, an inactive form Xq and 
an activated form X (often a phosphorylized version of 
Xq). The conversion of the two forms into each other is 
provided by an activating enzyme A (often a kinase), the 
deactivation by another enzyme D (often a phosphatase). 
In the activation process, the catalyst A first binds its 
substrate Xq- The resulting enzyme-substrate complex 
AXq may decay back into the original components. In 
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the case of a successful conversion, however, a product 
molecule X is released and the enzyme A is recovered for 
further use. The deactivation process is analogous. 

The CMC can be functionally decomposed into two 
enzymatic conversion processes. According to Michaelis- 
Menten kinetics (comp. Appendix VAl, the conversion 
rate is, in the linear regime, limited by the amount of 
available substrate. For very high substrate concentra- 
tion, however, the conversion rate approaches a max- 
imum value, determined only by the amount and effi- 
ciency of the enzyme (saturation regime). 

As demonstrated in a classical paper by Goldbeter and 
Koshland [Gol8l] , the combination of the two enzymatic 
conversion reactions can lead to interesting behavior if 
they operate within the saturated regime. In this case, 
the equilibrium ratio [A]/[Xo] as a function of the ratio 
of enzyme levels [A]/ [D] develops a sigmoidal shape with 
a sharp transition point (zero-order ultra-sensitivity) . In 
the context of biochemical signal networks, CMCs are for 
this reason understood as switches. 

The Goldbeter-Koshland theory is based on determin- 
istic (mass action) rate equations and thus disregards 
fluctuations entirely. Molecular reactions, however, in- 
evitably generate intrinsic noise, due to their discrete and 
stochastic nature. Even under so-called steady-state con- 
ditions, the momentary rates at which reactions proceed 
are fluctuating around the mean values described by mass 
action theory. The corresponding temporal fluctuations 
of molecule numbers are particularly important in living 
cells, where the average molecule numbers of many chem- 
ical species are low. For this reason, quantitative models 
of biochemical concentration fluctuations are developed 
for different types of reaction networks (see, for example, 
Refs. |Nac06l IWarOBl |Q5a02] ). 

Due to their ubiquity in living cells, CMCs are of par- 
ticular interest. A detailed theoretical investigation of 
the intrinsic fluctuations of CMCs, their robustness and 
tunability was provided by Levine et al. |Lev07] , who di- 
rectly solved the master equation for the probability dis- 
tribution of the number of activated signal molecules. 
The authors further consider the information transmis- 
sion properties of the system in the presence of the in- 
trinsic fluctuations, by applying a pulse-like increase of 
the kinase activity as an input signal. They find that the 
noisy CMC can transmit the signal reliably if tuned to 
an optimal parameter range. 

In this paper, we focus on the shape of the station- 
ary probability distributions produced by CMCs in var- 
ious parameter regimes. The reaction kinetics of this 
system is simulated using the exact Gillespie algorithm. 
This simulation yields directly the temporal concentra- 
tion fluctuations x(t) of the activated signaling molecule. 

We find an unexpected variety of distribution func- 
tions P(x), including Gaussian, exponential, flat, as well 
as power law distributions with a fractional and tunable 
exponent. The type of the emerging distribution func- 
tion depends on parameters such as the total amount of 
available enzyme and substrate molecules in their differ- 



ent forms and on reaction rate coefficients. We specu- 
late that living cells could switch between distinct sta- 
tistical distributions, on short time scales, by controlling 
the overall expression levels of these molecules. In many 
cases, moreover, the enzymes of a CMC are themselves 
activated and deactivated by another cycle. In this way, 
the effective conversion efficiency of an enzyme can be 
tuned over a wide range with only minimal changes of 
protein expression levels. This tremendous flexibility of 
CMCs with respect to their statistical properties suggests 
a more complex picture of cellular signal processing which 
is based on the active generation and precise shaping of 
concentration fluctuations of signaling molecules. 

In our paper we develop analytical approximations of 
the concentration fluctuations within CMCs, based on 
stochastic differential equations and explicit stationary 
solutions of the corresponding Fokker-Planck equations. 
The analytical results are in excellent agreement with the 
simulations and provide a quantitative understanding of 
the major statistical features. 



II. MODELS AND METHODS 
A. Model parameters and assumptions 

Let the reactions take place in a container of volume V, 
so that the concentration [S] of a substance corresponds 
to a molecule number s = [S]V. We also assume that the 
reactor is "well-stirred", i.e. diffusion of chemical species 
is infinitely fast and so spatial effects are disregarded. 

We study a CMC of the form 

61 c\ 
X + A . AX a >X + A 

Ul 

b 2 c 2 
X + D . DX >X + D. (1) 

The substrate X is converted into its activated form X 
by enzyme A. The corresponding deactivation is per- 
formed by enzyme D. We thus have to consider 6 tem- 
porally variable molecule numbers Xq, x, a, d, axg, dx, dy- 
namically coupled by 6 chemical reactions. Within each 
enzymatic conversion unit, the 3 reaction coefficients 
are denoted b (binding), u (unbinding) and c (conver- 
sion). Index 1 is used for the activation and index 2 
for the deactivation process. Additional parameters are 
the total amount of the substrate in its various forms, 
x t — xq + axQ + x + dx, as well as the total amounts of 
enzymes a f = a + axg and d t = d + dx. 



B. Analytical and numerical methods 

Analytical approximations for chemical reaction net- 
works can be obtained by deriving Langevin equations 
for the temporal changes of the molecule numbers. These 
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stochastic differential equations contain, besides a deter- 
ministic term that corresponds to the mass action change 
rates, a stochastic term that accounts for the fluctua- 
tions. To make use of the standard methods of stochas- 
tic differential calculus, the fluctuation term is approxi- 
mated by a Gaussian, white noise random process. This 
is a critical approximation, since the effective "strength" 
of the white noise process has to be chosen with care, 
in order to reflect the true process as faithfully as possi- 
ble. In the case of chemical Langevin equations, the true 
process consists of a series of delta-peaks, arriving with 
(inhomogcncous) Poisson statistics, ft is therefore possi- 
ble to derive the proper strength of the white noise pro- 
cess from the fundamental properties of Poisson statis- 
tics. This theory of chemical Langevin equations has 
been developed, for the general case, by Gillespie [GilOOJ. 
In this paper, we take a similar approach, suitable for 
our specific reaction network. 

In order to test our analytical approximations, we shall 
compare the results with a numerically exact Montc- 
Carlo-Simulation of the reaction dynamics by imple- 
menting the Gillespie algorithm |Gil77] . In this algo- 
rithm, the molecule numbers of each species are in- 
tegers which change abruptly due to elementary reac- 
tion events. Statistically, these elementary reactions 
are Poisson-processes with average event rates depend- 
ing on the momentary molecule numbers, according to 
the chemical rate equations. Therefore, the intrinsic 
stochastic fluctuations of the reactions are automatically 
included in a realistic way. 



C. Coarse graining of the enzymatic conversion 

We first focus on a single enzymatic conversion reac- 
tion, for example the activation process. Our goal is to 
describe it in a coarse grained approximation as a single 
functional unit with effective statistical properties. Two 
of these effective units will later be combined (as shown 
in Figjl]) to derive a stochastic differential equation for 
the fluctuating number x(t) of X- molecules. 

We assume for a moment that the number xq of sub- 
strate molecules X is constant (ideal reservoir). We are 
then interested in the average production rate R ac t(xo) 
of the activated protein X and in the temporal fluctua- 
tions AR act (xQ,t) of this rate. This, in turn, will enable 
us to write a stochastic rate equation of the production 
process in the form x = R a ct(xo) + AR act (xo 1 t). 

As for the average rates, we solve the mass action 
rate equations in the stationary state. This follows stan- 
dard Michaelis-Menten theory, but for completeness we 
include the derivation in Appendix |V A The result is 



Ract(xo) = U n 



•1'0 



Xq -\- k m 

with the maximum conversion velocity 



(2) 



v rn — ca t 

and the Michaelis constant 
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(3) 
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Enzymes are sometimes likened to nano machines, 
which convert their substrates in a predictable, goal- 
oriented process. Yet, many enzymes in biological sys- 
tems are working in a much more imperfect way: Once 
the enzyme has bound to its substrate, the enzyme- 
substrate-complex often dissociates back into the orig- 
inal two molecules. Each individual enzyme molecule 
will go through a series of futile binding-dissociation cy- 
cles, before it actually converts a substrate into the mod- 
ified form. In the chemical reaction equation 0, this is 
accounted for by the back reaction with rates Uj (with 
j = 1,2). The conversion efficiency of an enzyme can be 
quantified by the fraction of binding events that lead to 
a successful production and release of the modified sub- 
strate molecule. This fraction, in turn, depends on the 
relative magnitude of the rates Uj and Cj . We can define 
two limiting regimes: The case Uj » Cj corresponds to 
extremely inefficient enzymes. In the diagram of Fig. m, 
almost all activity of the reaction system will then take 
place within the shaded sub-modules. The flux in and 
out of these sub-modules is so weak that within the sub- 
modules a chemical equilibrium is established between 
the bound and dissociated enzyme-substrate-complexes. 
We therefore call this case the " pre-equilibrium" regime. 
The opposite case, Cj >> Uj corresponds to highly ef- 
ficient enzymes. In the diagram of Fig.Q, the system 
is running uni-directionally around the cycle, for most 
of the time. We therefore call this case the "sequential" 
regime. 

Independently from u and c, two other limiting regimes 
are connected with the amount of substrate xq relative to 
the Michaelis constant k m . The system is in the 'linear' 
regime for xq << k m and in the 'saturation' regime for 
>> k m . 

Next we model the fluctuations AR act (xa,t) of the pro- 
duction rate around the average value R a et(%o)- The 
statistical properties of these fluctuations are not obvi- 
ous, even if the substrate molecule number xo is arti- 
ficially held constant. As motivated in Appendix |V B 
we approximate the production process, in a coarse 
grained view, as a Poisson process with average event 
rate R a ct(xo)- Numerical simulations, shown below, con- 
firm that the probability distribution of the waiting time 
between successive X-production events is indeed ex- 
ponentially distributed with the expected characteristic 
time constant. We further approximate the above Pois- 
son process by white Gaussian noise with a proper pre- 



factor( Appendix VCl. As a result of the above coarse- 



graining procedure, we obtain 



X = Ract(xo) + V Ract(x a ) ■ £(*), 



(5) 
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where is normalized white Gaussian noise with 

(C(i)C(*')> = *(*-*')• 



D. Stochastic differential equation of a CMC 

We next combine the activation and deactivation pro- 
cesses. The molecule numbers x(t) and xo(t) are now 
both considered as variables. One obtains 



X = [ Ract(xo) - R dea (x)] 

\jRact(x Q ) ■ Ca(t) + \jRdea{x) ■ Q(t) 



(6) 



Note that the deactivation rates depend on x, not 
xq. To make further progress, we neglect the amount 
of substrates bound within enzyme-substrate complexes, 
so that xq — x t — x. Additionally, we make the sim- 
plifying assumption that the noise terms of the activa- 
tion and deactivation processes fluctuate statistically in- 
dependent from each other. We can then combine both 
terms, adding up the variances: 



Ract{x t -X) - Rdea{x)} 



Ract(x t ~X) + Rdea(x) 



CM- 



CO 



This has the general form of a stochastic differential 
equation with a multiplicative noise term [T] : 



Here, 



x = f(x) + g(x)-((t). 



,, , (x t -x) X 
f(x) = v a - r — ; v d - 



(x t ~ x) + k a x+ k d 



(8) 



(9) 



and 



g(x) 



(x t - x) 



'(x t -x) + k a ' d x + kd , 



(10) 



with obvious definitions of v a ,Vd,k a ,kd- In the follow- 
ing, we will extract statistical properties of this random 
process. Note that the Ito interpretation has to be used, 
whenever the true random process (that is to be approx- 
imated by Gaussian white noise) consists of a series of 
5— peaks, such as in our case of intrinsic, chemical noise 
|Ris84llKlm92] . 

We define a drift term, 



A(x) = f(x) 



B{x) = -g\x). 



(12) 



The time-dependent PDF P(x, t) of the fluctuating 
variable x(t) approximately satisfies the Fokker-Planck 
equation 



|p(M) = ~ [A(x)P(x,t)} + ^ [B(x)P(x,t)] 



(13) 



The stationary solution P(x) of this equation reads 



P(x) 



N 
Bjxj 



exp 



Ms) 
Ms) 



ds 



(14) 



Here, N is a normalization constant. 



E. The symmetric CMC 

With v a ,Vd,k a ,kd and Xt , there is obviously a large pa- 
rameter space to explore. In this paper, we shall restrict 
ourselves to just a few interesting cases. In a symmetric 
CMC, the activation and deactivation processes have the 
same parameters, i.e. v a = Vd — v and k a = kd = k. We 
then have 



fix) 



(x t 



and 



9 2 (x) 



(xt — x) + k x + k 



(x t - x) 



(15) 



(x* — x) + k x + k 



(16) 



Because the drift term A(s) and the diffusion term 
-B(s) are both proportional to v, it is clear that the max- 
imum production rate v will not affect the shape of the 
stationary PDF. Consequently, k and xt are the only im- 
portant parameters left. 



1. Linear Regime 

The limit of a large Michaclis constant, k >> xt, cor- 
responds to the linear regime of the two enzymatic con- 
version reactions. In this case, the terms x and (xt — x) 
can be neglected in Eqs. (15 1 and (16 1. This leaves us 
with 



f{x) = (vxt/k) - {2v/k)x 



(17) 



(11) and 



and a diffusion term 



9 2 {x) = (vxt/k). 



(18) 
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A straight forward calculation of the PDF yields a Gaus- 



sian, centered at x = with a variance <r\ = ^: 



ex e 



2(^-(x t /2))^ 



(19) 



The stochastic differential equation of a symmetric, lin- 
ear CMC corresponds to an Ornstein-Uhlenbeck process. 
Besides the Gaussian PDF, we therefore expect an expo- 
nentially decaying autocorrelation function: 



C xx (t) =< Ax(t)Ax(0) >= (^) 



,-{2v/k)r 



(20) 



The characteristic time constant is r r = 



2. Saturation Regime 

Next, we consider the opposite case of a small Michaelis 
constant, i.e. k « xt, corresponding to the saturation 
regime. We then have 



/(*) 



1 - 



vk 



for x » k 



and 



g 2 (x) = v 



x 

x + k 



2v for x >> k. 



(21) 



(22) 



The asymptotic drift and diffusion terms are A(x) — 
f , B{x) = v, and A(s)/B(s) = |. Therefore, 



Ms) 



ds = k- log(x/x min ), 



(23) 



and 



P(x) cx e ^os(x/ Xmm ) K ( X / Xm . n )k_ 



(24) 

Hence, we expect an increasing power-law tail for the 
asymptotic PDF in the saturation regime of the symmet- 
ric CMC. The exponent of the power-law can be frac- 
tional and is equal to the dimensionless Michaelis con- 
stant (Eq|l|. The above analytical approximations will 
break down when x approaches the limits or i t . 



F. The asymmetric CMC 

We now allow the activation parameters k a and v a to 
differ from the corresponding deactivation parameters kd 
and Vd- Under saturation conditions (x t » k a , xt » 
kd) and in the limit of large x one obtains f(x) — > (v a — 
Vd) and g 2 {x) — > (v a + Vd), so that 

A(s) 



A =2- 



(25) 



B(s) v a + v d 

This results in a stationary PDF with an exponential tail: 

P(x) cx e Aa: . (26) 

The decay constant A is positive for v a > Vd and negative 
for v a <v d . 



III. RESULTS 
A. Validation of Poisson statistics 

We first investigate the statistics of the enzymatic ac- 
tivation process, with artificially fixed number xq of sub- 
strate molecules. For this purpose, we perform direct 
Monte-Carlo simulations in different parameter regimes. 
All rates and times are presented in dimensionless num- 
bers. 
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FIG. 2: (Color online) Monte-Carlo simulation of enzymatic 
conversion: Molecule numbers of the enzyme-substrate com- 
plex (solid) and of the activated product (dashed) in the case 
of only one enzyme molecule. Parameters: 6 = u = c = 1.0, 
et = 1. The vertical arrows denote a conversion (c), binding 
(b) and unbinding (u) process. At is the time interval between 
two successive conversion events. 

The stochastic time evolution of the enzymatic acti- 
vation process is characterized by abrupt changes of the 
various molecule numbers by integer amounts (Fig. |2|. A 
single enzyme molecule sometimes undergoes binding (6) 
and unbinding (u) without conversion (c) to a product 
molecule. The time interval At between two successive 
conversion events is fluctuating around the inverse of the 
average production rate. 

Since the production process involves a sequence of el- 
ementary reaction steps, the distribution function P(At) 
of this waiting time is not expected to be exponential for 
an individual enzyme molecule. However, the superpo- 
sition of many such multi-step processes running inde- 
pendently from each other can closely mimic a Poisson 
process (Fig|3|. 



B. Monte-Carlo Simulation of the CMC 

Next we discuss the statistical properties of CMCs in 
selected parameter regimes, as obtained by Monte-Carlo 
simulation of the reaction dynamics. We shall mainly 
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FIG. 3: (Color online) Waiting time distributions: Simu- 
lated PDF of the time intervals between successive conver- 
sion events. Parameters: b = u = c = 1.0, xo = 1 = const.. 
Case (a): Only one enzyme molecule. Case (b): 10 indepen- 
dent enzyme molecules. The inset shows the same data in a 
semi- logarithmic plot. 



focus on CMCs with symmetric parameters for the ac- 
tivation and deactivation process. The total number of 
substrate molecules x* was 200 in all cases. Our analytic 
theory was based on the assumption that the amount 
of substrate bound in complexes is small compared to 
x t . We have therefore chosen a small number of enzyme 
molecules, a t = d t =10. The (rounded) parameters for 
all following simulations are listed in Tab.(|l]}. We have 
also included a saturation parameter (SP), defined as 
SP = x t /k m , and an equilibrium parameter (EP), de- 
fined as EP = c/u. For instance, SP > l,EP > 1, 
would indicate that the system is in the saturated, pre- 
cquilibrium regime. 
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ci = 1.5 
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ci = 1.25 
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1 


0.1 


1 


ci = 1.125 
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0.1 


1 




1.1 


180 


0.1 



TABLE I: . Parameter space explored in Monte-Carlo simu- 
lations. 



1. Symmetric CMC in the linear and weakly saturated 
regimes 

In the linear regime, we expect for the substrate X 
a Gaussian distribution, peaked at x = x*/2 and with 
variance xt/4. The autocorrelation of the random vari- 
able x(t) should decay exponentially with time constant 
t c = k/2v. The agreement of the Monte-Carlo results 
with this analytic theory is excellent (see Fig.Q). In the 
weakly saturated regime, we find a decrease of the aver- 
age molecule number and a considerable broadening of 
the distribution, while the shape of the PDF remains ap- 
proximately Gaussian. The distributions do not change 
dramatically when the parameter regime is changed from 
sequential to pre-equilibrium conditions, as long as the 
ratio of enzyme to substrate molecules is small (see foot- 
note 0). 




Number of X molecules 



FIG. 4: (Color online) X- distributions in the linear regime 
(b,c) and in the weakly saturated regime (d,e). The solid line 
(a) is the analytical solution to the linear case. Inset: Normal- 
ized auto- correlation function for linear case (symbols) with 
analytical solution (solid line). Parameters see Tab. 0). 

The Monte-Carlo simulations also yield the distribu- 
tions of the enzyme molecule number E = a = d (see 
Fig.(|5|). The effect of sequential or pre-equilibrium con- 
ditions is almost invisible for the particular parameters 
chosen. 



2. Symmetric CMC in the saturation regime 

Next, we turn to CMCs operating within the satu- 
ration regime, which corresponds to the hypersensitive, 
'switch-like' mode of the cycle. In the simulations, k m 
was indirectly changed via the conversion rate c. While 
small conversion rates result in a Gaussian PDF, the dis- 
tributions become extremely asymmetric as the system 
runs into the saturation regime (Fig. [6]). The double- 
logarithmic plot reveals a power law wing at the 'left' 
side of the peak. The positive exponent of the power law 
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FIG. 5: (Color online) Semilog arithmic plots of the simu- 
lated enzyme distributions in the linear regime (a,b) and in 
the weakly saturated regime (c,d). Parameters see Tab. Q). 



tail is fractional in the general case. It is determined by 
the Michaelis constant, as expected from the analytical 
theory above. For a very small Michaelis constant, one 
obtains an almost flat distribution, which can cover sev- 
eral decades of concentration. Of course, the PDF has 
sharp cutoffs at the maximum particle number x = x t 
and close to x — (not shown). 

This remarkable result demonstrates that the notion 
of deterministic biomolccular networks, with well-defined 
average levels of concentration and negligibly small Gaus- 
sian fluctuations, dramatically fails in certain parame- 
ter ranges. Concentration fluctuations with a power law 
wing are scale- free, and therefore arbitrarily large devi- 
ations from the average value occur with non-negligible 
probability. 




10 10 
Number of X molecules 

FIG. 6: (Color online) Pouier law tails in the strongly satura- 
tion regime: Double-logarithmic plots of the X-distribution. 
In cases (a)-(e) the conversion rate c has been gradually in- 
creased. Parameters see Tab. (III). 



3. Asymmetric CMC in the saturation regime 

From a systems biology point of view, an interesting 
question is the sensitivity of the CMC with respect to 
its parameter values. In particular, we investigated the 
effects of tuning the system slightly away from the com- 
pletely symmetric parameter settings considered so far. 
The most dramatic effects are expected for a CMC in the 
hypersensitive saturation regime. 

For this purpose, we start again with the parameters of 
the symmetric saturated CMC, which produced a PDF 
with a power law tail of slope 1.1 (compare Fig. |6je)). 
Now, however, we fine-tune the conversion rate c\ of the 
activation reaction, while leaving the corresponding pa- 
rameter C2 at its former value 1. 

As expected, if c 2 < c l5 the PDF of X is peaked 
around a small average concentration, while X has a 
high average concentration (Fig. [7]). The average con- 
centrations are drastically different even for rather sim- 
ilar c-parameters, due to the hypersensitive response of 
the saturated CMC. We find PDFs with exponential tails 
for all cases, except in a very narrow range around per- 
fect parametric symmetry. This is in agreement with the 
analytical theory presented in section |IIF| In the nar- 
row symmetrical regime, the two PDFs collapse to one. 
They are mirror-symmetric with respect to the average 
molecule number in this case. 



S 10" 

'8 

B 

i 



-a 

x> 
o 

Oh 



10 



-2 



10 



"(a) 




(a) 


; \(b) 




(b)A 




■' v> 

' :\\(c) 
1 1 * 

ii * 


(d) 


(cyyAl 




\A > 

A 1 

f ' 1 \ 
t ' ' * 




/ /v 1 

/ Ik 





50 100 150 200 
Number of molecules 



FIG. 7: (Color online) Collapse of exponential distributions at 
the critical point of parametric symmetry: Distributions of Xq 
(dashed lines) and X (solid lines). In cases (a)-(d) the con- 
version rate of the activating reaction only has been gradually 
increased. The left wing of (d) corresponds to Figplc) when 
plotted double-logarithmically. Parameters see Tab. Jill. 

This behavior somewhat resembles critical phenom- 
ena in physics, where fluctuations of arbitrary size occur 
when a control parameter is precisely tuned to a critical 
value. 

In biological systems, it would be extremely improb- 
able to find a CMC where all the microscopic parame- 
ters of the activation and deactivation reaction are pre- 



cisely identical. However, equations (15 1 and (16 1 show 
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that effective dynamical symmetry can be achieved un- 
der the much weaker conditions v a sa Vd and k a w kd- In 
terms of the microscopic parameters, this translates into 
dot « c 2 d t and » aga. 

In order to demonstrate that the dynamics is only con- 
trolled by the conversion velocities, the Michaelis con- 
stants and the total amount of substrate Xt, we have 
performed another Monte-Carlo simulation for a CMC 
with x t = 200 and the microscopically non-symmetric 
parameters at — 50, b\ = 1, u\ = 0.1, ci = 1 for the activa- 
tion reaction and d t — 10, 62 = 10, — 6, C2 = 5 for the 
deactivation reaction. These parameters are neverthe- 
less symmetrical on the coarse-grained level of k and v. 
The simulation results indeed show a power law behavior, 
thus confirming the analytical prediction (see Fig. pi). 
Note that the total amount of enzymes a t and d t can 
be easily varied in a living cell, for example by chang- 
ing the expression levels or the activity of the enzymes. 
This offers a way to tune the CMC through the critical 
point. If, for instance, we detune a t away from the crit- 
ical value af rlt%> = 50 by ±10 percent, we find that one 
of the distributions P(X) and P(Xq) is loosing its power 
law behavior. Yet, the (respective) complementary form 
of substrate still shows a very steep power law tail under 
these conditions of disturbed symmetry. 



IV. DISCUSSION AND OUTLOOK 

The statistical properties of concentration fluctuations 
produced by CMCs reveal an extremely rich behavior. A 
variety of qualitatively different probability distributions 
has been found for the molecule numbers of the acti- 
vated substrate, depending on the parameter settings. A 
particularly remarkable result for symmetric CMCs op- 
erated in the saturation regime was the emergence of an 
extremely broad PDF with a power law tail. These fluc- 
tuations are driven by purely intrinsic noise, originating 
from the stochastic arrival times of the molecular reac- 
tion events. 

We note that in biological systems there are additional, 
extrinsic sources of noise as well. For example, we have 
considered the total number of enzyme molecules, at and 
d t , as being strictly constant in this paper. In biological 
systems, the enzymes are themselves subject to produc- 
tion and consumption processes and will therefore un- 
dergo concentration fluctuations. When these enzymes 
serve a CMC in the saturation regime, the steady state 
activation level x /xq of the substrate will depend hyper- 
scnsitivcly on the momentary ratio of enzyme concen- 
trations a t /d t . Small (and sufficiently slow) fluctuations 
of the enzyme concentrations will therefore be amplified, 
leading to an additional, extrinsic broadening of the PDF 
of x(t). 

At first glance, it seems that such extreme concentra- 
tion fluctuations would compromise the function of bio- 
chemical networks [ThaOll ITha02j . However, recent re- 
ports have suggested that large biochemical fluctuations 




Number of molecules 



FIG. 8: (Color online) Tuning the CMC through the critical 
point by changing the enzyme concentration a t . The simulated 
CMC has asymmetric rate constants, but becomes symmetric 
with respect to the effective coarse-grained parameters k m and 
v m for at = 50. Parts (a)-(c) correspond to a t = 45, 50 and 55. 
Shown are double-logarithmic distributions of the activated 
(solid lines) and deactivated (dashed lines) substrates. For 
parameters see text. 

can also be beneficial for organisms, ranging from bac- 
teria to humans |Aus061 lRao02l IHasOOl IPauOOj . In a 
recent review article [Los08 , Losick and Desplan have 
summarized a number of studies showing that certain 
cells choose one or another pathway of differentiation 
stochastically, without regard to environment or history. 
Another example of stochastic signal processing is pro- 
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vided by the well-understood bacterial chemotaxis net- 
work. The flagellar motor of the bacterium is normally 
rotating in the counterclockwise (CCW) direction, but 
shows stochastic intervals of clockwise (CW) rotation. 
This gives rise to distinct phases of straight swimming 
motion of the bacterium, separated by random tumbling 
phases. Cell-membrane receptors detect the concentra- 
tion of attractant molecules in the surrounding medium 
of the bacterium. Over several intermediate steps, the ac- 
tivation level of the receptors affects the distribution of 
CCW interval length and, thereby, the run length distri- 
bution of the bacterium's random walk in the medium. 
A statistical analysis of the CCW intervals revealed a 
power-law distribution jKor04| . which has been related 
to molecular noise in the reaction network |Tu05| . Inter- 
estingly, such random walks with power law-distributed 
run lengths (Levy-flights) are known to generate trajec- 
tories which are the optimum strategies to search effi- 
ciently for randomly located objects |Vis931 . This exam- 
ple shows how the shaping of molecular noise and the 
modulation of the noise parameters in response to envi- 
ronmental stimuli can be used by cells for complex tasks, 
such as foraging behavior. 

We note that similar ideas of stochastic signal process- 
ing have recently emerged in the field of neuro science 
|Gro09] . In the new concept of 'reservoir computing', a 
network of (randomly) connected neurons generates a so- 
called transient state dynamics, where the trajectory of 
the system state is temporally fluctuating between vari- 
ous unstable attractors. This autonomously active 'reser- 
voir' network is only weakly coupled to the 'input' and 
'output' units. As simulations have shown, the mapping 
of low dimensional input signals onto the high dimen- 
sional state space of the reservoir network can be advan- 
tageous for the signal processing. 

Finally, in this report we have discussed the station- 
ary behavior of a single CMC in which the total number 
of molecules is fixed. In living cells, however, multiple 
CMCs are connected in linear and branched signaling 
networks. Moreover, the total number of molecules fluc- 
tuates as new proteins are expressed or old proteins are 
recycled. If already a single CMC under stationary con- 
ditions gives rise to such highly complex, bizarre and 
non-deterministic behavior as described in this report, 
we argue that concentration fluctuations in living cells 
are even less predictable by classical mass action theory. 



V. APPENDIX 

A. Average production rates 

We consider an enzymatic conversion reaction of the 
general form: 



b c 
X + E t - Y >Z + E 



Using mass action rate theory, we obtain for the temporal 
change of the concentration y(t) of the enzyme-substrate 
complex: 



y = b x e — u y — cy 



(28) 



We make the simplifying approximations that x(t) is 
held constant. After a certain relaxation time, the system 
will reach a steady state, in which also y(t)=const. The 
condition y = then leads to 



x e- 



(29) 



The expression = is defined as the inverse 
Michaelis constant, so that y = Since the enzyme 
can either be free or bound in the complex, e t — e + y, 

k m ■ 



one obtains y - 



Solving for y yields 



V = e t 



3C ~\~ k<n 



(30) 



For the quantity of interest, the steady state generation 
rate z = cy of the product, we finally obtain 



(c e t ) 



(31) 



(27) 



B. Enzymatic conversion as an effective Poisson 
process 

In general, an individual A-enzyme molecule can un- 
dergo a series of binding/unbinding events with the (non- 
exhaustible) substrate Xq, before the substrate is fi- 
nally converted into a new A'-molcculc. Therefore, even 
though each elementary reaction step, i.e. binding, un- 
binding and conversion, is a Poisson process, the same is 
not true for the multi-step production process. [3] 

However, many individual A-cnzyme molecules, dis- 
persed throughout the volume of the container, are si- 
multaneously active, with independent temporal statis- 
tics. Our numerical simulations show that the super- 
position of many independent non-Poisson processes can 
resemble an effective Poisson process very closely. As ex- 
pected, the characteristic time constant of this effective 
Poisson process is given by the inverse of the average 
total production rate R ac t{xo). 

In our CMC system, the substrate molecule number xq, 
and therefore R a ct(xo), are not constant. The resulting 
Poisson process is therefore not stationary but has a time- 
varying rate. 

We conclude that in systems with many independent 
enzyme molecules, the overall conversion process can be 
approximated by an inhomogeneous Poisson process. 
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C. Poisson process as white Gaussian noise 



Assume now a Poisson process with constant average 
event rate k = R act (xo). We express the temporal change 
of the number x(t) of product molecules in the form 



x = k + Ak(t). 



(32) 



For later convenience, we want to approximate the fluc- 
tuation term by Gaussian white noise, 



(AfcO)Afc(i')) = k5{t-t'). 



Dividing this equation by k leads to 



\ Vk Vk 



(37) 



(38) 



(Ak(t)Ak(t')) = T5(t-t'). 



(33) 



What is the proper choice for the pre- factor T, so that 
the major statistical properties of a Poisson process are 
consistently reproduced ? 

To answer this question, we consider the number n(T) 
of X-moleculcs which are produced during an interval of 
length T: 

n(T)= [ x(t)dt = kT + [ Ak(t)dt = n + An. (34) 
Jo Jo 

In the ensemble average, a Poisson process must fulfill 



((An) 2 ) 



(35) 



We now define a new stochastic process by 



C(*) = 



Ak(t) 



(39) 



It is also normally distributed, but shows the desired 
property of ^-autocorrelation with unit strength:. 



<C(i)C(*')> =*(*-*')• 



(40) 



We conclude that a proper description of a Poisson 
process by a stochastic differential equation should have 
the form 



or 



Y 

Ak(t)dt 



kT. 



(36) 



The left side of the above equation can be reduced 
to TT. Using Eq.(33l, we therefore obtain T = k, and 
therefore 



x = k+Vk((t). 
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[1] Note that stochastic differential equations of the general 
form x = f(x) + g(x) ■ £(i) are extremely rich in behavior 
and can produce random fluctuations with arbitrary PDF 
and ACF, as shown i n Re f. [Pri99] and |Pri00] 

[2] Note that in section HE we have neglected the amount 
of substrate which is bound in complexes. In order to 
refine the theory, let us define a new dynamic variable 
a = x + dx (The complementary variable f3 — xo + axo 
is unnecessary, since (3 — x t — a). This variable a defines 
the macro-state of the system in our coarse-grained view. 
It is changed only by activation or deactivation processes. 
On the other hand, binding and unbinding processes only 
affect the micro-state of the system. The latter is defined 
by the numbers dx and axo, each of which can vary be- 
tween and the respective number of enzyme molecules. 
Thus, each macro state a can be sub-divided into several 
micro states (dx,axo). The fluctuations of our variable of 
interest, x(t), are determined by changes of the macro- 
and of the micro-state. In the pre-equilibrium regime, for 
each momentary macro-state a, we expect that equilib- 
rium distributions P eq (dx\a) (and P eq (axo\/3)) of micro- 
states are building up. The probability of having x acti- 
vated substrate molecules is under such conditions given 
by P(x) = E Q >, P(<*) P«i{dx = a-x\a). 

[3] For a simple example, consider a sequence of one bind- 
ing and one conversion step. The PDF of each elementary 
Poisson step is exponential. The PDF of the sequence is 
a convolution of two exponential functions, i.e. a Gamma 
distribution with shape parameter k = 2. 



